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Abstract. We investigate the quantum state of burst atoms seen in the recent Rb-85 
experiments at JILA. We show that the presence of a resonance scattering state can 
lead to a pairing instability generating an outflow of atoms with energy comparable 
to that observed. A resonance effective field theory is used to study this dynamical 
process in an inhomogeneous system with spherical symmetry. 
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1. Introduction 

The ability to dynamically modify the nature of the microscopic interactions in a Bose- 
Einstein condensate — an ability virtually unique to the field of dilute gases — opens the 
way to the exploration of a range of fundamental phenomena. A striking example of 
this is the Bosenova experiment carried out in the Wieman group at JILA pQ which 
explored the mechanical collapse instability arising from an attractive interaction. This 
collapse resulted in an unanticipated burst of atoms, the nature of which is a subject 
of current debate. In this article we suggest a possible mechanism for the formation of 
these bursts by application of an effective quantum field theory which includes explicitly 
the resonance scattering physics. 

The Bosenova experiment conducted by the JILA group consisted of the following 
elements. A conventional stable Bose-Einstein condensate was created in equilibrium. 
The group then utilized a Feshbach resonance to abruptly switch the interactions to be 
attractive inducing an implosion. One might have predicted that the rapid increase in 
density would simply lead to a rapid loss of atoms, primarily through inelastic three- 
body collisions. In contrast, what was observed was the formation of an energetic burst 
of atoms emerging from the implosion. Although the energy of these atoms was much 
larger than that of the condensate, the energy was insignificant when compared to 
the molecular binding energy which characterizes the energy released in a three-body 
collision. In the end, what remained was a remnant condensate which appeared distorted 
and was believed to be in a highly excited collective state. 

One theoretical method which has been extensively explored to explain this 
behavior has been the inclusion of a decay term into the Gross-Pitaevskii equation as a 
way to account for the atom loss El HI E] • Aside from its physical application to the 
Bosenova problem, the inclusion of three-body loss as a phenomenological mechanism 
represents an important mathematical problem, since the nonlinear Schrodinger 
equation allows for a class of self-similar solutions in the unstable regime. The local 
collapses predicted in this framework can generate an outflow even within this zero 
temperature theory. However, there are a number of aspects which one should consider 
when applying the extended Gross-Pitaevskii equation to account for the observations 
made in the JILA experiment. 

The first problematic issue is the potential breakdown of the principle of attenuation 
of correlations. This principle is essential in any quantum or classical kinetic theory 
as it allows multiparticle correlations to be factorized. This assumption is especially 
evident in the derivation of the Gross-Pitaevskii equation where all explicit multiparticle 
correlations are dropped. However, as shown in Figure ^ even a simple classical model 
may exhibit clustering when mechanically unstable which appears to invalidate the 
assumption of an attenuation of correlations. Furthermore, there is also considerable 
evidence for this instability toward pair formation in the mechanically unstable quantum 
theory 

A second difficulty with motivating the Gross-Pitaevskii approach is that by this 
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Figure 1. A collapsing classical system with point-like objects interacting via a 
Gaussian attractive potential. The top graph shows the initial configuration. The lower 
graph displays the evolved state demonstrating clustering of the particles (illustrated 
as yellow regions). The red lines mark the particle trajectories. 



method one describes the interactions as energy independent through a single parameter, 
the scattering length, which is determined from the s-wave scattering phase shift at 
zero scattering energy. Near a Feshbach resonance, the proximity of a bound state in 
a closed potential to the zero of the scattering continuum can lead to a strong energy 
dependence of the scattering. Exactly on resonance, the s-wave scattering length passes 
through infinity and, in this situation, the Gross-Pitaevskii equation is undefined. 

These two fundamental difficulties with the Gross-Pitaevskii approach led us to 
reconsider this problem. We were motivated by the fact that the same experimental 
group at JILA recently performed a complementary experiment [Zj which provided key 
insights into the Bosenova system. What was remarkable in these new experiments 
was that, even with a large positive scattering length in which the interactions were 
repulsive, a burst of atoms and a remnant condensate were observed. Furthermore, in 
the large positive scattering length case, simple effective field theories which included an 
explicit description of the Feshbach resonance physics were able to provide an accurate 
quantitative comparison with the data [HI UD] • The theory showed the burst to arise 
from the complex dynamics of the atom condensate coupled to a coherent field of exotic 
molecular dimers of a remarkable physical size and near the threshold binding energy. 

In this paper, we draw the connections between the two JILA experiments. We 
pose and resolve the question as to whether the burst of atoms in the Bosenova 
collapse could arise in a similar way as in the Ramsey fringe experiment — from the 
formation of a coherent molecular superfluid. This hypothesis is tested by applying an 
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Figure 2. A Feshbach resonance results when a closed channel potential possesses a 
bound state in proximity to the scattering energy in an open channel potential. The 
detuning of the bound state from the edge of the collision continuum is denoted by vq. 

effective field theory for resonance superfluidity to the collapse. For fermions, the case 
of resonance superfluidity in an inhomogeneous system has been treated in the local 
density approximation, using essentially the uniform solution at each point in space 
[TT] . For the collapse of a Bose-Einstein condensate, as we wish to treat here, the local 
density approximation is not valid, and the calculation must be performed on a truly 
inhomogeneous system. This represents the first time that the resonance superfluidity 
theory has been applied to a system of this type. 

2. Effective Field Theory 

In the Feshbach resonance illustrated in Figure the properties of the collision of 
two ground state atoms is controlled through their resonant coupling to a bound state 
in a closed channel Born-Oppenheimer potential. By adjusting an external magnetic 
field, the scattering length can be tuned to have any value. This field dependence of 
the scattering length is characterized by the detuning uq and obeys a dispersive profile 
given by a(u Q ) = Ob g (l — k/(2v )), with k the resonance width and ab g the background 
scattering length. In fact, all the scattering properties of a Feshbach resonance system 
are completely characterized by just three parameters U = 4:7rh 2 a hg /m, g = \f kUq, 
and uq. Physically, Uq represents the energy shift per unit density on the single particle 
eigenvalues due to the background scattering processes, while go, which has dimensions 
of energy per square-root density, represents the coupling of the Feshbach resonance 
between the open and closed channel potentials. 

We now proceed to construct a low order many-body theory which includes this 
resonance physics. The Hamiltonian for a dilute gas of scalar bosons with binary 
interactions is given in complete generality by 

H = Jd 3 x ^ a {x)H a {x)^ a {x) + J d 3 xd 3 x' ^ a (x)^ a (x')V(x, x')^ a {x')^ a {x) (1) 
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where H a (x) is the single particle Hamiltonian, V(x,x') is the binary interaction 
potential, and ip a ( x ) is a bosonic scalar field operator. In cold quantum gases, where the 
atoms collide at very low energy, we are only interested in the behavior of the scattering 
about a small energy range above zero. There exist many potentials which replicate the 
low energy scattering behavior of the true potential; therefore, it is convenient to carry 
out the calculation with the simplest one, the most convenient choice being to take the 
interaction potential as a delta-function pseudopotential when possible. 

For a Feshbach resonance this choice of pseudopotential is generally not available 
since the energy dependence of the scattering implies that a minimal treatment must 
at least contain a spread of wave-numbers which is equivalent to the requirement of 
a nonlocal potential. Since the solution of a nonlocal field theory is inconvenient, we 
take an alternative but equivalent approach. We include into the theory an auxiliary 
molecular field operator ip m ( x ) which obeys Bose statistics and describes the collision 
between atoms in terms of two elementary components: the background collisions 
between atoms in the absence of the resonance interactions and the conversion of atom 
pairs into molecular states. This allows us to construct a local field theory with the 
property that when the auxiliary field is integrated out, an effective Hamiltonian of the 
form given in Eq. ((!} is recovered with a potential V(x, x') = V(\x— x'\) which generates 
the form of the two-body T- matrix predicted by Feshbach resonance theory [12] . The 
local Hamiltonian which generates this scattering behavior: 



has the intuitive structure of resonant atom-molecule coupling. Here V a ^ m are the 
external potentials and /x a m are the chemical potentials. The subscripts a, m represent 
the atomic and molecular contributions, respectively. The Feshbach resonance is 
controlled by the magnetic field which is incorporated into the theory by the detuning 
^ = Hm — Z^a between the atomic and molecular fields. The Hamiltonian in equation (J2J) 
contains the three parameters U, g, and v which account for the complete scattering 
properties of the Feshbach resonance. It is important to keep in mind that they are 
distinct from the bare parameters Uq, go, and uq introduced above. In order for the local 
Hamiltonian given in equation @ to be applicable, one must introduce into the field 
theory a renormalized set of parameters each containing a momentum cutoff associated 
with a maximum wavenumber K. This need not be physical in origin but should 
exceed the momentum range of occupied quantum states. The relationships between the 
renormalized and bare parameters are given by: U = TUo, g = Tg , and v = vo + aggo/2, 
where a = mK / '(2n 2 h 2 ) and T = (1 — aUo)^ 1 [T2]- All the results presented here have 
been shown to be independent of the momentum cutoff in the theory. 




(2) 
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We define the condensates in terms of the mean-fields of the operators 4> a (x) — 
(ijj a (x)) and <f) m (x) = (tp m (x)) along with the fluctuations about the atomic field 
Xa(x) = i]) a {x) — 4> a (x). Note that, in principle, there is also a term which involves 
the fluctuations about 4> m (x). Assuming the occupation of 4> m (x) to be small (less than 
2% in the simulations we present), we drop higher order terms arising from fluctuations 
about this mean-field which do not give a significant correction to our results. As 
discussed in reference [TJ] , we derive four equations: two corresponding to a Schrodinger 
evolution of the mean fields 

lh ^M_ = (_|lv2 + K(^)-/ia + f/[|0a( : r)| 2 + 2G w (a ; , : r)])0 a ( : r) 

+ [UG A (x,x)+g(f> m (x))(f>* a {x) 
lh ^M. = (-^-Vl + Vjx) - ii m )4>m{v) + + G A (x, x)) (3) 

and two corresponding to the Louiville space evolution of the normal density 
Gn(x,x') = (xl(x')xa(x)) and of the anomalous density Ga(x,x') = (xa(x')xa(x)) 

%h d -^ = Y,Q-QY). (4) 
The density matrix and self-energy matrix are defined respectively as 

[ (xl(x')Xa(x)) (Xa(x')Xa(x)} \ 

[ (xl(x')xt(x)) (Xa(x')xi(x)) ) 



S{x,x') 



H(x x')- ( H ^ X>) ^ X>) \ (5) 

H x,x)- y _ A ^ x/) _ H * {XjX/) )■ W 

The convenience of choosing a microscopic model in which the potential couplings are of 
contact form is now evident since the elements of the self-energy matrix S are diagonal 
in x and x' with non-zero elements 

h 2 

H(x, x)= - —W 2 X + V a (x) -n a + 2U[\<p a (x)\ 2 + G N (x, x)] 

A(x, x) = U[<t> 2 a (x) + G A (x, x)} + g<f> m (x). (6) 

Since the normal density and anomalous pairing field are both six-dimensional objects, 
it is very difficult to solve these equations in an arbitrary geometry. For this reason we 
consider the case of greatest symmetry consisting of a spherical trap. Here we can reduce 
the problem to one of only three dimensions, which is still nontrivial to treat, according 
to the following procedure. To begin with, it is convenient to write the elements of the 
single particle density matrix in center of mass R = (x + x')/2 and relative r = x — x' 
coordinates. The normal density then takes on a familiar structure corresponding to 
the Wigner distribution 

G N (R, k)= j d 3 r (xl(R - r/2) Xa (R + r/2))e~ ik - r 

d 3 rG N {R,r)e~ ik - r (7) 



which in the high-temperature limit will map on to the particle distribution function 
f(R, k) for a classical gas. Correspondingly, the anomalous density can also be written 
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Figure 3. Illustration of the spherically symmetric geometry used, as defined by the 
center of mass vector R, relative momentum vector fc, and the angle 9 between them. 



in Fourier space as 

Ga(R, k)= J d 3 r ( Xa (R ~ r/2) Xa (R + r/2)) e - iW 

= J d 3 rG A {R,r)e- ik - r . (8) 

In this geometry, the angular dependence of R is irrelevant, and the cylindrical 
symmetry about R allows the wavevector k to be represented by its length and the 
one remaining angle as illustrated in Figure El This allows us to represent the density 
distributions in three-dimensions as 

G{R,k) = G{R,k,6) (9) 

where G corresponds to either the normal (Gn) or anomalous (Ga) density. It is now 
straightforward to rewrite equations Q and (J3J) in this coordinate system. It is worth 
pointing out the simple structure of the kinetic energy contributions to equation (jlj) 
which for the GV and Ga components take the corresponding forms respectively: 

(Vl - Vl) G N (x, x') = 2 (Vfl ■ V r ) Gn(R, r) (10) 

1, 



{yl + vl) g a (x, x') 
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One may now take the Fourier transform with respect to r as indicated by equations U\ 
and |Hl replacing V r — > ik. The gradient operator V# can be expressed in any 
representation, but it is most convenient to use spherical polar coordinates aligned with 
the k direction vector 

Vfl = fl4 + ( jI4- + 4 3J_» (12) 



where tp is the azimuthal angle about k (which will eventually drop out in our chosen 
symmetry), and R, 8, and (p are the spherical unit vectors in the R, 6, and ip directions, 
respectively. Noting that R - k = kcos9, 9 ■ k — —ksin9, and (p ■ k = 0, we arrive at 
the following expression for the differential operator in equation ITOl 



V R -k 



' d sin 9 d " 
COS OR ~ ~R~~d6 



(13) 
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Furthermore, the spherical Laplacian for a system with no azimuthal dependence, as 
required in equation ^2 is given by 

^9 1 d / „ 9 d \ 1 d f . „ d \ , J , 

v « = jpag( B as) + jpsswH*)- (14) 

In practice, we expand the 9 dependence of Gn and Ga in terms of the orthogonal 
Legendre polynomials, and the angular derivatives are then easily implemented via the 
usual recursion relations. 

3. Results and Analysis 

As an initial test, we expect the resonance theory to give a similar prediction to the 
Gross-Pitaevskii equation in the initial phase of the collapse when the quantum depletion 
is small. Figure 0] shows a direct comparison between the Gross-Pitaevskii approach 
and the resonance theory. The same initial conditions were used for all our simulations; 
1000 rubidium-85 atoms in the ground state of a 10 Hz harmonic trap. For all the 
images we present, the results of the three-dimensional calculation, in our spherical 
geometry, are illustrated as a two-dimensional slice through the trap center. In the 
Gross-Pitaevskii solution we used a scattering length of —200 a where ao is the Bohr 
radius. For comparison, the Feshbach resonance theory uses a positive background 
scattering length of 50 ao and a resonance width and detuning respectively of 15 kHz 
and 2.8 kHz. These parameters give the same effective scattering length as the one used 
in the Gross-Pitaevskii evolution, but nowhere in the resonance theory does the effective 
scattering length appear explicitly. As is evident, there is no noticeable discrepancy 
between the two approaches over this short timescale. Eventually we expect these 
theories to diverge significantly as the density increases and the coupling between the 
atomic and molecular degrees of freedom becomes stronger. However, at this stage the 
agreement is a demonstration that our renormalized theory correctly allows us to tune 
the interactions in an inhomogeneous situation. 

We now proceed to a more complex situation in which the timescales for the atom- 
molecule coupling and the collapse dynamics are more compatible. From a numerical 
point of view, it becomes convenient to increase the resonance width to 1.5 MHz and 
the detuning to 14 kHz so that the effect of the atom-molecule coupling will appear in 
the first stage of the collapse. This allows us to form a complete picture of the dynamics 
involving the atomic collapse and the simultaneous coupling to a coherent molecular 
field. The numerical calculation is shown as a movie in Fig. |3] for both the condensed 
and noncondensed components. One sees the formation of a significant fraction of 
noncondensed atoms-a feature not described within the Gross-Pitaevskii framework. 
During a time evolution of 0.8 ms the condensate fraction falls to approximately 80% of 
its initial value, while the noncondensate fraction reaches a peak at around 20%. The 
amplitude of the scalar field m remains below the 2% level at all times. 

To better illustrate the behavior of the atoms during the collapse we present the 
flow of the different distributions involved. The condensate velocity field is shown in 
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Figure 4. A direct comparison of the collapse between the Gross-Pitaevskii (left) 
and the resonance approach (right) within the regime of applicability of the Gross- 
Pitaevskii equation. Each horizontal pair is at the same time step with time increasing 
from top to bottom. As expected we observe no appreciable difference between the 
two methods. 

Figure |3 It exhibits similar characteristics to those predicted by the Gross-Pitaevskii 
theory, which without loss, predicts that the condensed atoms will always accelerate 
toward the trap center. In contrast, the velocity field of the atoms outside of the 
condensate is radially outward. The production of this component is quite interesting 
because it is this same component which in the theory of the Ramsey fringe experiment 
was quantitatively determined to give rise to the burst. 

Obviously an important quantity to calculate for these expanding noncondensed 
atoms is the effective temperature, or energy per particle, since this quantity is observed 
experimentally. This is illustrated in Figure where superimposed on an illustration 
of the density is a colormap of the temperature. The hottest atoms generated in the 
center of the cloud are of comparable energy scale to that seen in the experiment, being 
on the order of 100 nK. 

4. Conclusion 

We emphasize that the work presented here is a model calculation to illustrate the 
feasibility for the burst to be generated through atom-molecule coupling. However, 
there are a number of important distinctions with the experimental situation which 
would have to be accounted for before making a direct comparison. These simulations 
contain no inelastic three-body loss and particle number is absolutely conserved. In 
reality three-body loss may be important to the experiment, but we suggest with this 
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Figure 5. The simulation of the collapse in the resonance theory showing the time 
evolution of the condensed fraction <j) a {x) (left) and noncondensed fraction G_n(x,x) 
(right). Each horizontal pair is taken at the same instant of time, with time increasing 
from top to bottom. It is evident that non condensate atoms are produced during 
the collapse dynamics, forming rings which propagate from the center of the cloud 
outward. 

work that three-body loss is not the only mechanism for producing a noncondensed 
burst during the collapse. 

It should be emphasized that if our hypothesis for the burst generation is correct, the 
noncondensate atoms that are produced by this mechanism are not simply generated in 
a thermal component, but are instead generated in a fundamentally intriguing quantum 
state. The process of dissociation of molecules into atom pairs produces macroscopic 
correlations reminiscent of a squeezed vacuum state in quantum optics. This means 
that every atom in the burst with momentum k would have an associated partner with 
momentum —k. In principle, the correlations could be directly observed in experiments 
through coincidence measurements providing clear evidence as to whether this is the 
dominant mechanism for the burst generation in the Bosenova. 
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Figure 6. The velocity fields for the condensate component (f) a (x) (left) and the 
noncondensed component Gn(x,x) (right) midway through the simulation (0.4 ms). 
The color contours indicate the densities and the velocity fields are represented in 
direction and strength by the arrows. This clearly shows, that in the resonance theory 
as the condensate collapses inward, the noncondensate atoms that are generated flow 
outward. 
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Figure 7. Density distribution of the noncondensate atoms near the end of the 
simulation (0.8 ms) on which we have superimposed the energy per particle as a 
colormap. The range of energies, of order 100 nK, is consistent with the characteristic 
scale of the burst particle energies in the Bosenova experiment. Note that hot atoms 
are generated in the center of the cloud during the atom-molecule oscillations since 
this is where the atom-molecule coupling is strongest (the coupling strength varies as 
the square root of the density). As the hot particles radiate outward a ring can be 
observed. 
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